﻿#include "canny.h"

int Wsrgb2ycc[12] = { 58,149,13,4096,
                    -31,-81,112,32768,
                    112,-103,-9,32768 };
int Wsycc2rgb[12] = { 298,0,431,-59995,
                    298,-48,-167,22789,
                    298,550,0,-75227 };
int* wsp;

inline int GetRGB2YCC_Y(int r, int g, int b) {
    wsp = &Wsrgb2ycc[0];
    return clip((wsp[0] * r + wsp[1] * g + wsp[2] * b + wsp[3]) >> 8, 0, 255);
}

inline int GetRGB2YCC_Cb(int r, int g, int b) {
    wsp = &Wsrgb2ycc[4];
    return clip((wsp[0] * r + wsp[1] * g + wsp[2] * b + wsp[3]) >> 8, 0, 255);
}

inline int GetRGB2YCC_Cr(int r, int g, int b) {
    wsp = &Wsrgb2ycc[8];
    return clip((wsp[0] * r + wsp[1] * g + wsp[2] * b + wsp[3]) >> 8, 0, 255);
}

inline int GetYCC2RGB_R(int y, int cb, int cr) {
    wsp = &Wsycc2rgb[0];
    return clip((wsp[0] * y + wsp[1] * cb + wsp[2] * cr + wsp[3]) >> 8, 0, 255);
}

inline int GetYCC2RGB_G(int y, int cb, int cr) {
    wsp = &Wsycc2rgb[4];
    return clip((wsp[0] * y + wsp[1] * cb + wsp[2] * cr + wsp[3]) >> 8, 0, 255);
}

inline int GetYCC2RGB_B(int y, int cb, int cr) {
    wsp = &Wsycc2rgb[8];
    return clip((wsp[0] * y + wsp[1] * cb + wsp[2] * cr + wsp[3]) >> 8, 0, 255);
}

inline char RGB2YCC(int r, int g, int b, int& y, int& cb, int& cr) {
    y = GetRGB2YCC_Y(r, g, b);
    cb = GetRGB2YCC_Cb(r, g, b);
    cr = GetRGB2YCC_Cr(r, g, b);
    return OK;
}

inline char YCC2RGB(int y, int cb, int cr, int& r, int& g, int& b) {
    r = GetYCC2RGB_R(y, cb, cr);
    g = GetYCC2RGB_G(y, cb, cr);
    b = GetYCC2RGB_B(y, cb, cr);
    return OK;
}




/*
* convolution
*/
char Conv2D(matrix inmat, matrix outmat, matrix gaustable) {
    if (inmat.w <= 0 || inmat.data == NULL || outmat.w <= 0 || outmat.data == NULL ||
        gaustable.w <= 0 || gaustable.d == NULL) {
        printf("<GetEdge_Gradient>ERROR! inputs is empty!\n");
        return FALSE;
    }
    int r, c, ir, ic, chan, index, sum = 0;
    if (gaustable.sum <= 0) {                   //get sum of table
        for (r = 0; r < gaustable.h; r++)
            for (c = 0; c < gaustable.w; c++)
                sum += ((int*)gaustable.d)[r * gaustable.w + c];
        gaustable.sum = sum;
    }
    for (r = 0; r < inmat.h; r++) {
        for (c = 0; c < inmat.w; c++) {
            for (chan = 0; chan < inmat.channel; chan++) {
                int t = 0;
                if (r < (gaustable.h >> 1) || r >= inmat.h - (gaustable.h >> 1) ||
                    c < (gaustable.w >> 1) || c >= inmat.w - (gaustable.w >> 1)) {
                    outmat.data[(r * outmat.w + c)* outmat.channel + chan] = (uchar)t;
                    continue;
                }
                for (ir = 0; ir < gaustable.h; ir++) {
                    for (ic = 0; ic < gaustable.w; ic++) {
                        index = ((r + ir - (gaustable.h >> 1)) * inmat.w + c + ic - (gaustable.w >> 1)) * inmat.channel + chan;
                        t += inmat.data[index] * ((int*)gaustable.d)[ic * gaustable.w + ir];
                    }
                }
                t /= gaustable.sum;                
                if (t < 0 || t>255) {
                    printf("<GetEdge_Gradient>ERROR! outof range!\n");
                    return FALSE;
                }
                outmat.data[(r * outmat.w + c)* outmat.channel + chan] = (uchar)t;
            }

        }
    }
    return OK;
}

/*
* Get Grayimg
*/
uchar *GetGrayimg(matrix img) {
    uchar* p = NULL, * head = NULL, * q = NULL;
    int r, c;
    head = p = (uchar*)malloc(img.w * img.h);
    if (head == NULL) {
        printf("<GetGrayimg>ERROR! not enough memory!\n");
        return head;
    }
    q = img.data;
    for(r=0;r<img.h;r++)
        for (c = 0; c < img.w; c++) {
            *p = int(q[0] * 0.299 + q[1] * 0.587 + q[2] * 0.114); //Y = 0.299 * R + 0.587 * G + 0.114 * B
            q += 3;
            p++;
        }
    return head;
}

/**************************************
* like GetAngle2 
* 
**************************************/
inline int GetAngle3(int dx, int dy) {
    //dy = abs(dy);
    float dy_dx = (int)abs(dy * 1000 / (((double)dx) + 1e-5));

    if (abs(dx) < 1)
        dy_dx = 9999999;
    float tan225 = 0.414213562373, tan675;
    tan675 = 1.0 / tan225;
    if(dy_dx<=tan225&&dy_dx>=-tan225)
        return DIRECTIONMARK::DRIGHT;
    if (dy_dx <= tan675 && dy_dx >= tan225)
        return  DIRECTIONMARK::DRIGHTDOWN;
    if (dy_dx > tan675 || dy_dx < -tan675)
        return DIRECTIONMARK::DDOWN;
    if (dy_dx <= -tan225 && dy_dx >= -tan675)
        return DIRECTIONMARK::DLEFTDOWN;
    return DIRECTIONMARK::DRIGHT;
}

/**************************************
* cal angle
* take the positive half axis of x as the polar axis
* accuracy=90.0/L   //精度
* 'x' exlarged 1000 times
**************************************/
int GetIndextheta(int x, int* arr, int l) {
    int i;
    for (i = 0; i < l; i++) {
        if (x < *arr)
            return i;
        arr++;
    }
    return l;
}
/**************************************
* Get direction2
*     \     |     /
*      1    2   3 
*        \  |  /  
* -----0----------0right--------
*        /  |  \  
*3leftdown 2down 1rightdown
*     /     |      \
* 'direction mark ' see enum of DIRECTIONMARK
**************************************/
inline int GetAngle2(int dx, int dy,cannypara cfg) {
    int tanv = (int)abs(dy * 1000 / (((double)dx)+1e-7));    
    if (abs(dx) < 1)
        tanv = 9999999;
    float theta = 0;
    int lv = 0;
    if (dx >= 0) {
        if (dy >= 0) {      //first quadrant
            theta = (float)(GetIndextheta(tanv, cfg.at, cfg.atl) * 90.0 / cfg.atl);
        }
        else {              //fourth quadrant
            theta = 360-(float)(GetIndextheta(tanv, cfg.at, cfg.atl) * 90.0 / cfg.atl);
        }
    }    else {
        if (dy >= 0) {      //second quadrant
            theta = 180-(float)(GetIndextheta(tanv, cfg.at, cfg.atl) * 90.0 / cfg.atl);
        }
        else {              //third quadrant
            theta = (float)(GetIndextheta(tanv, cfg.at, cfg.atl) * 90.0 / cfg.atl) + 180;
        }
    }
    lv = ((int)((theta > 180 ? theta-180 : theta)+22.5))/45;     //cal direction
    if (lv > 4 || lv < 0) {
        printf("<GetAngle2>error!\n");
        exit(-1);
    }
    lv = (lv == 4 ? 0 : lv);
    //debug....
    //printf("dx:%d\tdy:%d->tan:%d\t", dx, dy, tanv);
    //printf("theta:%.2f\tlv:%d\n", theta,lv);

    return lv;
}
/****************************************
* Get direction 
*        \         |         /
*           \   1  |  2    /
*        0     \   |    /  3
* -----------------------------------
*   leftdown  /    |    \   right
*           /      |       \
*        / down    |rightdown\
***************************************/
inline int GetAngle(int dx,int dy) {
    //dy = abs(dy);
    if (dx >= 0) {
        if (dy >= 0) {
            if (dx > dy)
                return DIRECTIONMARK::DRIGHT;
            else
                return DIRECTIONMARK::DRIGHTDOWN;
        }
         else
            if (dx > -dy)
                return DIRECTIONMARK::DLEFTDOWN;
            else
                return DIRECTIONMARK::DDOWN;

    }
    if (dy > -dx)
        return DIRECTIONMARK::DDOWN;
    return DIRECTIONMARK::DLEFTDOWN;
}
inline int GetEdgeLv(int x, int ths[2]) {
    if (x < ths[0])
        return 0;    
    if (x < ths[1])
        return 1;
    return 2;
}

/*
* test some function
* 
* 
*/
char TestCANNY(cannypara cfg) {
    int t=0;
    t += GetAngle2(2, 0, cfg);       //lv---0
    t += GetAngle2(20, 1, cfg);       //lv---0
    t += GetAngle2(20, -1, cfg);       //lv---0
    t += GetAngle2(-100, 1, cfg);       //lv---0
    t += GetAngle2(-100, -1, cfg);       //lv---0
    t += GetAngle2(-100, 0, cfg);       //lv---0
    t += GetAngle2(2, 2, cfg);       //lv---1
    t += GetAngle2(20, 18, cfg);       //lv---1
    t += GetAngle2(19, 20, cfg);       //lv---1
    t += GetAngle2(-2, -2, cfg);       //lv---1
    t += GetAngle2(-20, -18, cfg);       //lv---1
    t += GetAngle2(-19, -20, cfg);       //lv---1
    t += GetAngle2(0, 3, cfg);       //lv---2
    t += GetAngle2(1, 13, cfg);       //lv---2
    t += GetAngle2(-1, 23, cfg);       //lv---2
    t += GetAngle2(0, -3, cfg);       //lv---2
    t += GetAngle2(2, -43, cfg);       //lv---2
    t += GetAngle2(-3, 53, cfg);       //lv---2
    t += GetAngle2(-13, 13, cfg);       //lv---3
    t += GetAngle2(-14, 13, cfg);       //lv---3
    t += GetAngle2(-17, 18, cfg);       //lv---3
    t += GetAngle2(13, -13, cfg);       //lv---3
    t += GetAngle2(14, -13, cfg);       //lv---3
    t += GetAngle2(17, -18, cfg);       //lv---3
    if (t == 36) {
        printf("test 'GetAngle2' success!\n");
        return OK;
    }
    else {
        printf("test 'GetAngle2' failed!\n");
        return FALSE;
    }
}
/*
* Get edge map and angle map
* input：
* gtable-gaussian noise filter
* cfg-sobel parameter   (including angle table)
* 
*/
char GetEdgeGaussSobel(matrix grayimg,matrix edgemap,matrix anglemap, matrix gtable,cannypara cfg) {
    int r, c, i, j;
    int tx[3][3] = { {-1,0,1},{-2,0,2},{-1,0,1} };
    int ty[3][3] = { {-1,-2,-1},{0,0,0},{1,2,1} };
    int x, y, t,a;
    
    matrix img;
    memcpy(&img, &grayimg, sizeof(matrix));
    img.data = NULL;
    img.data = (uchar*)malloc(img.w * img.h * img.channel);
    if (img.data == NULL) {
        printf("<GetEdgemapSobel>ERROR! out of memory!\n");
        return FALSE;
    }
    memcpy(img.data, grayimg.data, img.w * img.h * img.channel);
    Conv2D(img, grayimg, gtable);
    for (r = 0; r < grayimg.h; r++) {
        for (c = 0; c < grayimg.w; c++) {
                x = y = 0;
                //if (r == 1727 && c == 922) {        //debug...
                //    x = 0;
                //}
                if (!(r<1 || r>grayimg.h - 2 || c<1 || c>grayimg.w - 2))
                    for (i = 0; i < 3; i++) {
                        for (j = 0; j < 3; j++) {
                            t = grayimg.data[(r + i - 1) * grayimg.w + c + j - 1];
                            y += ty[i][j] * t;
                            x += tx[i][j] * t;
                            //printf("%dx%d=%d\t", ((int*)tablex.d)[index],t,x);
                        }
                        // std::cout << std::endl;
                    }
                //printf("sum:%d,t:%d", sum, t);
                x = x/4;
                y = y/4;
                t = (int)sqrt((x * x) + (y * y));
                // printf("res:%d\n", t);
                if (t < 0 || t>255) {
                    printf("<GetEdgemapSobel>warning! out of range t=%d\n", t);
                }
                //t = GetEdgeLv(clip(t, 0, 255), cfg.dt)*127;
                a = GetAngle2(x, y,cfg) * 85;
                edgemap.data[r * edgemap.w + c] = (uchar)(clip(t, 0, 255));
                anglemap.data[r * anglemap.w + c] = (uchar)(a);            
        }
    }

    free(img.data);
    return OK;
}

/*
* non-maximum suppression
* 
* flag :0---neighbor
* flag :1---single direction
*/
char NonMaxSupp(matrix dstmat,matrix srcedgemap, matrix anglemap,char flag) {
    if (dstmat.w <= 0 || dstmat.data == NULL || srcedgemap.w <= 0 || srcedgemap.data == NULL ||
        anglemap.w <= 0 || anglemap.data == NULL) {
        printf("<NonMaxSupp>ERROR! inputs is empty!\n");
        return FALSE;
    }
    int r, c, ir, ic, d1,d2,e1,e2, size = 3,index,book=0;
    // compare 3x3 neighbor
    if (flag == 0) {
        for (r = 0; r < srcedgemap.h; r++) {
            for (c = 0; c < srcedgemap.w; c++) {
                if (r < 1 || r >= srcedgemap.h - 1 || c < 1 || c >= srcedgemap.w - 1) {
                    continue;
                }
                if (srcedgemap.data[r * srcedgemap.w + c] == 0)
                    continue;
                d1 = anglemap.data[r * anglemap.w + c];
                e1 = srcedgemap.data[r * srcedgemap.w + c];
                book = 0;
                for (ir = 0; ir < size; ir++) {
                    if (book > 0)
                        break;
                    for (ic = 0; ic < size; ic++) {
                        index = (r + ir - (size >> 1)) * srcedgemap.w + c + ic - (size >> 1);
                        d2 = anglemap.data[index];
                        if (d2 == d1) {         //same diretion
                            e2 = srcedgemap.data[index];
                            if (e2 > e1) {
                                dstmat.data[r * srcedgemap.w + c] = 0;
                                book++;
                                break;
                            }
                        }
                    }
                }
            }
        }
    }
    else {
        int e0, d0;        //single direction compare 1x3 or 3x1 or... neighbor
        for (r = 0; r < srcedgemap.h; r++) {
            for (c = 0; c < srcedgemap.w; c++) {
                if (srcedgemap.data[r * srcedgemap.w + c] == 0)
                    continue;
                if (r < 1 || r >= srcedgemap.h - 1 || c < 1 || c >= srcedgemap.w - 1) {
                    continue;
                }
                d0 = anglemap.data[r * anglemap.w + c];
                e0 = srcedgemap.data[r * srcedgemap.w + c];
                if (d0 < 85) {          //diretion 0
                    e1 = srcedgemap.data[r * srcedgemap.w + c - 1];
                    e2 = srcedgemap.data[r * srcedgemap.w + c + 1];
                }
                else if (d0 < 170) {    //diretion 1
                    e1 = srcedgemap.data[(r - 1) * srcedgemap.w + c - 1];
                    e2 = srcedgemap.data[(r + 1) * srcedgemap.w + c + 1];
                }
                else if (d0 < 255) {    //diretion 2
                    e1 = srcedgemap.data[(r - 1) * srcedgemap.w + c];
                    e2 = srcedgemap.data[(r + 1) * srcedgemap.w + c];
                }
                else {  //diretion 3
                    e1 = srcedgemap.data[(r - 1) * srcedgemap.w + c + 1];
                    e2 = srcedgemap.data[(r + 1) * srcedgemap.w + c - 1];
                }
                if (e1 > e0 || e2 > e0) {
                    dstmat.data[r * srcedgemap.w + c] = 0;
                }
            }
        }
    }
    return OK;
}

char DoubleThre(matrix dstmat, matrix srcedge, cannypara cfg) {
    int r, c,lv;
    for (r = 0; r < srcedge.h; r++) {
        for (c = 0; c < srcedge.w; c++) {            
            lv = GetEdgeLv(srcedge.data[r * srcedge.w + c], cfg.dt);
            dstmat.data[r * srcedge.w + c] = lv*127;
        }
    }
    return OK;
}


char EdgeTrick(matrix dstmat, matrix srcmat,int size=3,char flag=0) {
    int r, c, ir, ic, x, y, book;
    for (r = 0; r < srcmat.h; r++) {
        for (c = 0; c < srcmat.w; c++) {
            if (srcmat.data[r * srcmat.w + c] != WEAKEDGE*127)
                continue;
            book = NONEDGE;
            if(c==1172&&r==1847)
                book = NONEDGE;
            for (ir = 0; ir < size; ir++) {
                if (book == 1)
                    break;
                for (ic = 0; ic < size; ic++) {
                    x = c + ic - (size >> 1);
                    y = r + ir - (size >> 1);
                    if (x < 0 || x >= srcmat.w || y < 0 || y >= srcmat.h)
                        continue;
                    if (srcmat.data[y * srcmat.w + x] == STRONGEEDGE * 127) {
                        book = STRONGEEDGE;
                        break;
                    }
                }
            } 
            dstmat.data[r * srcmat.w + c] = book * 127;
        }
    }

    return OK;
}

/*
* my canny
* 
*/
char mycanny(matrix inmat, matrix outmat, matrix gaustable, cannypara cfg){
    if (inmat.w == 0 || inmat.h == 0 || inmat.data == NULL) {   //check data...
        printf("<mycanny> input is empty!\n");
        return FALSE;
    }
    if (outmat.w == 0 || outmat.h == 0||outmat.data==NULL) {
        printf("<mycanny> save is empty!\n");
        return FALSE;
    }
    if (gaustable.w == 0 || gaustable.h == 0 || gaustable.d == NULL) {
        printf("<mycanny> input gaustable is empty!\n");
        return FALSE;
    }
    int r, c, sum = 0;

    for (r = 0; r < gaustable.h; r++)                    //get sum of table
        for (c = 0; c < gaustable.w; c++)
            sum += ((int*)gaustable.d)[r * gaustable.w + c];
    gaustable.sum = sum;
    
    matrix grayimg;
    grayimg.h = inmat.h; grayimg.w = inmat.w; grayimg.channel = 1;
    grayimg.data = NULL;
    if (inmat.channel == 3) {                   //Get gray image
        grayimg.data = GetGrayimg(inmat);
        if (grayimg.data == NULL) {
            printf("<mycanny>ERROR! can not get grayimg!\n");
            return FALSE;
        }
    }    else {
        grayimg.data = (uchar*)malloc(grayimg.w * grayimg.h);
        if (grayimg.data == NULL) {
            printf("<mycanny>ERROR! not enough memory!\n");
            return FALSE;
        }
        memcpy(grayimg.data, inmat.data, inmat.w * inmat.h);
    }

    matrix edgemap; matrix anglemap,nonmaxmat,dtmat,trackmat;
    edgemap.w = anglemap.w = grayimg.w;
    edgemap.h = anglemap.h = grayimg.h;
    edgemap.channel = anglemap.channel = 1;
    edgemap.data = anglemap.data = NULL;
    edgemap.data = (uchar*)malloc(edgemap.w * edgemap.h);
    anglemap.data = (uchar*)malloc(anglemap.w * anglemap.h);
    memcpy(&nonmaxmat, &edgemap, sizeof(matrix));
    memcpy(&trackmat, &edgemap, sizeof(matrix));
    memcpy(&dtmat, &edgemap, sizeof(matrix));
    nonmaxmat.data = (uchar*)malloc(edgemap.w * edgemap.h);
    trackmat.data = (uchar*)malloc(trackmat.w * trackmat.h);
    dtmat.data = (uchar*)malloc(dtmat.w * dtmat.h);
    if (edgemap.data == NULL || anglemap.data == NULL|| nonmaxmat.data==NULL||
        trackmat.data==NULL|| dtmat.data==NULL) {
        printf("<mycanny>not enough memory!\n");
        return FALSE;
    }
    memset(edgemap.data, 0, edgemap.w* edgemap.h);
    memset(anglemap.data, 0, anglemap.w* anglemap.h);

    TestCANNY(cfg);

    //main process
    if (OK != GetEdgeGaussSobel(grayimg, edgemap, anglemap, gaustable, cfg))
        return FALSE;
    memcpy(nonmaxmat.data, edgemap.data, edgemap.w* edgemap.h);
    if (OK != NonMaxSupp(nonmaxmat, edgemap, anglemap, 0))
        return FALSE;

    DoubleThre(dtmat, nonmaxmat, cfg);

    memcpy(trackmat.data, dtmat.data, trackmat.w * trackmat.h);
    if (OK != EdgeTrick(trackmat, dtmat))
        return FALSE;
    memcpy(outmat.data, inmat.data, outmat.w * outmat.h* outmat.channel);
    BilateraMask(outmat, inmat, gaustable, trackmat);

    //debug...
    writeMat(grayimg, "./mark/gaussiagrayimg.ppm", PPM_P6_OUT);
    writeMat(edgemap, "./mark/edgeimg.ppm", PPM_P6_OUT);
    writeMat(anglemap, "./mark/angleimg.ppm", PPM_P6_OUT);
    writeMat(nonmaxmat, "./mark/nonmaxmat.ppm", PPM_P6_OUT);
    writeMat(dtmat, "./mark/doublethreshold.ppm", PPM_P6_OUT);
    writeMat(trackmat, "./mark/trackmat.ppm", PPM_P6_OUT);

    free(edgemap.data);
    free(anglemap.data);
    free(grayimg.data);
    free(nonmaxmat.data);
    free(trackmat.data);
    free(dtmat.data);
    return OK;
}

/*
* Do 'Bilatera' after 'DoubleThre'
*  
* 
*/
char NonMaxDoubleThre(matrix inmat, matrix outmat, matrix gaustable, cannypara cfg) {
    if (inmat.w == 0 || inmat.h == 0 || inmat.data == NULL) {   //check data...
        printf("<NonMaxDoubleThre> input is empty!\n");
        return FALSE;
    }
    if (outmat.w == 0 || outmat.h == 0 || outmat.data == NULL) {
        printf("<NonMaxDoubleThre> save is empty!\n");
        return FALSE;
    }
    if (gaustable.w == 0 || gaustable.h == 0 || gaustable.d == NULL) {
        printf("<NonMaxDoubleThre> input gaustable is empty!\n");
        return FALSE;
    }
    int r, c, sum = 0;
    for (r = 0; r < gaustable.h; r++)                    //get sum of table
        for (c = 0; c < gaustable.w; c++)
            sum += ((int*)gaustable.d)[r * gaustable.w + c];
    gaustable.sum = sum;
    matrix grayimg;
    grayimg.h = inmat.h; grayimg.w = inmat.w; grayimg.channel = 1;
    grayimg.data = NULL;
    if (inmat.channel == 3) {                   //Get gray image
        grayimg.data = GetGrayimg(inmat);
        if (grayimg.data == NULL) {
            printf("<NonMaxDoubleThre>ERROR! can not get grayimg!\n");
            return FALSE;
        }
    }
    else {
        grayimg.data = (uchar*)malloc(grayimg.w * grayimg.h);
        if (grayimg.data == NULL) {
            printf("<NonMaxDoubleThre>ERROR! not enough memory!\n");
            return FALSE;
        }
        memcpy(grayimg.data, inmat.data, inmat.w * inmat.h);
    }

    matrix edgemap; matrix anglemap, nonmaxmat, dtmat,masktextbg, maximg;
    edgemap.w = anglemap.w = grayimg.w;
    edgemap.h = anglemap.h = grayimg.h;
    edgemap.channel = anglemap.channel = 1;
    edgemap.data = anglemap.data = NULL;
    edgemap.data = (uchar*)malloc(edgemap.w * edgemap.h);
    anglemap.data = (uchar*)malloc(anglemap.w * anglemap.h);
    memcpy(&nonmaxmat, &edgemap, sizeof(matrix));
    memcpy(&dtmat, &edgemap, sizeof(matrix));
    memcpy(&masktextbg, &edgemap, sizeof(matrix));
    memcpy(&maximg, &inmat, sizeof(matrix));
    nonmaxmat.data = (uchar*)malloc(edgemap.w * edgemap.h);
    dtmat.data = (uchar*)malloc(dtmat.w * dtmat.h);
    masktextbg.data = (uchar*)malloc(dtmat.w * dtmat.h);
    maximg.data = (uchar*)malloc(maximg.w * maximg.h* maximg.channel);
    if (edgemap.data == NULL || anglemap.data == NULL || nonmaxmat.data == NULL ||
        dtmat.data == NULL|| masktextbg.data==NULL) {
        printf("<mycanny>not enough memory!\n");
        return FALSE;
    }
    memset(edgemap.data, 0, edgemap.w * edgemap.h);
    memset(anglemap.data, 0, anglemap.w * anglemap.h);

    TestCANNY(cfg);

    writeMat(grayimg, "grayimg.ppm", PPM_P6_OUT);
    //main process
    if (OK != GetEdgeGaussSobel(grayimg, edgemap, anglemap, gaustable, cfg))
        return FALSE;
    memcpy(nonmaxmat.data, edgemap.data, edgemap.w * edgemap.h);
    if (OK != NonMaxSupp(nonmaxmat, edgemap, anglemap, 1))
        return FALSE;
    DoubleThre(dtmat, nonmaxmat, cfg);

    //debug...
    writeMat(grayimg, "mark\\gaussiagrayimg.ppm", PPM_P6_OUT);
    writeMat(edgemap, "mark\\edgeimg.ppm", PPM_P6_OUT);
    writeMat(anglemap, "mark\\angleimg.ppm", PPM_P6_OUT);
    writeMat(nonmaxmat, "mark\\nonmaxmat.ppm", PPM_P6_OUT);
    writeMat(dtmat, "mark\\dtmat.ppm", PPM_P6_OUT);

    double wts[3] = { 1.0,0.5,0 };
    GetMaskTextBg(inmat, masktextbg, dtmat);
    BilateraMasklvs(outmat, inmat, gaustable, dtmat, masktextbg,wts,3);

    writeMat(outmat, "res.ppm", PPM_P6_OUT);

    addimgs(inmat, dtmat, maximg);
    writeMat(maximg, "mark/maximg.ppm", PPM_P6_OUT);
    writeMat(masktextbg, "mark/masktextbg.ppm", PPM_P6_OUT);

    free(edgemap.data);
    free(anglemap.data);
    free(grayimg.data);
    free(nonmaxmat.data);
    free(dtmat.data);

    return OK;
}

char GetMaskTextBg(matrix inmat, matrix outmat,matrix edgemap) {
    if (inmat.w != outmat.w || inmat.data == NULL || outmat.data == NULL) {
        printf("input is error!\n");
        return FALSE;
    }
    int r, c, chan, y, cb, cr, ir, ic, xi, yi, y0, cb0, cr0;
    int kernelsize = 5;
    int thre_max=10,deltay=30,ymax=150;

    uchar* p=inmat.data,*pout=outmat.data;
    int count, sumv,index,colors,edgeflag;
    cv::Mat img=cv::Mat(edgemap.h,edgemap.w,CV_8UC1);
    img.data = edgemap.data;
    
    for (r = 0; r < inmat.h; r++) {
        for (c = 0; c < inmat.w; c++) {
           /* if (edgemap.data[r * edgemap.w + c] == 0) {
                *pout = 0;
                pout++;
                p += 3;
                continue;
            }*/
            if (c == 1190 && r == 1945)
                edgeflag = 0;
            edgeflag = 0;
            for (ir = 0; ir < 3; ir++) {
                for (ic = 0; ic < 3; ic++) {
                    yi = r + ir - (3 >> 1);
                    xi = c + ic - (3 >> 1);
                    if (xi < 0 || xi >= edgemap.w || yi < 0 || yi >= edgemap.h) {
                        continue;
                    }
                    if (edgemap.data[yi * edgemap.w + xi] != 0) {
                        edgeflag = 1;
                        break;
                    }
                }
                if (edgeflag != 0)
                    break;
            }
            if (edgeflag == 0) {
                *pout = 0;
                p += 3;
                pout++;
                continue;
            }
            count = sumv = 0;
            RGB2YCC(p[0], p[1], p[2], y0, cb0, cr0);
            if (abs(cb0-128) + abs(cr0-128) < thre_max&&y0< ymax) {  //black text ,,Maximum color component
                for (ir = 0; ir < kernelsize; ir++) {
                    for (ic = 0; ic < kernelsize; ic++) {
                        yi = r + ir - (kernelsize >> 1);
                        xi = c + ic - (kernelsize >> 1);
                        if (xi < 0 || xi >= inmat.w || yi < 0 || yi >= inmat.h) {
                            continue;
                        }
                        else {
                            index = (yi * inmat.w + xi)* 3;
                            RGB2YCC(inmat.data[index], inmat.data[index+1], inmat.data[index+2], y, cb, cr);
                            colors = abs(cb - 128) + abs(cr - 128);
                            if (abs(y - y0) < deltay&&colors< thre_max) {      //like source point(black text)
                                continue;
                            }
                            count++;
                            sumv += colors;
                        }
                    }
                }
            }
            if (count == 0) {
                *pout = 0;
            }
            else {
                if ((sumv / count) < thre_max) {        //black edge on white background
                    int t;
                    t = 255 - int(sumv / count);
                    *pout = (uchar)clip(t, 0, 255)>50?255:0;
                    //printf("out value:%d\n", *pout);
                }
                else {
                    *pout = 0;
                }
                //*pout = (uchar)clip(int(sumv / count) <<4, 0, 255);
            }
            pout++;
            p += 3;
        }
    }
    writeMat(outmat, "mark/BkEdge_WBg.ppm", PPM_P6_OUT);

    return OK;
}




